Summary of the manuscript “On the possibilities of predicting geomagnetic 
secular variation with geodynamo modeling” to be submitted to Geophysical 
Research Letter (GRL) 

It has long been known that the Earth possesses an internal magnetic field (i.e. the field vanishes at infinite 
distance from the Earth). This geomagnetic field is believed to be generated and maintained by convective 
flow in the Earth’s liquid outer core (geodynamo). In die past decade, several numerical models have been 
developed to model the geodynamo, including our MoSST (Modular, Scalable, Self-consistent and Three- 
dimensional) core dynamics model. These models can successfully explain qualitatively the geodynamo 
process in the outer core, e.g. a dominantly dipolar geomagnetic field at the surface (with the polarity 
almost parallel to die geometric polarity of the Earth), westward-drift of geomagnetic field lines and 
reversals of geomagnetic polarity. 

However, no attempt has been given to quantitative applications of the numerical models on geomagnetic 
field, partly due to large differences between the parameters used in numerical simulations and those 
appropriate for the Earth’s core, and partly due to physical approximations adopted in numerical modeling. 
But such studies are important for geodynamo and geomagnetic research: surface geomagnetic 
observations can be used to constrain numerical geodynamo models, and numerical models can be applied 
to forecast geomagnetic secular variation observable on and near the Earth’s surface. 

This research article reports our first ever effort on applying our MoSST core dynamics model and the 
observed surface geomagnetic field to predict geomagnetic secular variation. The surface geomagnetic 
field is obtained via the comprehensive field model operated here in GSFC. As the first attempt, we focus 
on examining how numerical dynamo solutions are affected by geomagnetic observation. For this purpose, 
a pure numerical geodynamo solution is selected to be assimilated with the surface geomagnetic field in 
1940. The assimilation is simple: the observed field is inserted into the dynamo solution. The new solution 
is then used as an initial state for numerical simulation. The simulated solutions are then used to compare 
with the observed surface geomagnetic field in the subsequent years. 

Our findings are very encouraging. While there is no correlation between the field from pure dynamo 
simulation and the field from observation, the field from the new solutions with data assimilation, in 
particular the large-scale (or low degree) field evolves closely with the observed field over time. As the 
result, the ass imila ted field can capture large-scale features, such as south Atlantic anomaly observed at the 
surface of the Earth for the period from 1940 to 1990. However, there are still discrepancies between the 
small-scale (high degree) field. In addition, the assimilated solutions diverge (without further assimilation 
constraint) above approximately 60-year periods. 

Our research results suggest that it is possible to assimilate numerical results with surface observations to 
predict future changes in geomagnetic field. They also suggest that further research is necessary on better 
understanding the statistical properties of numerical dynamo solutions, in particular the error development 
in time, and the stability of the numerical solutions under arbitrary perturbations. These are necessary for 
future development/implementation of better assimilation technologies. 

This research is supported by NASA Solid Earth and Natural Hazard Program (SENH) and by NSF 
Geophysics and Mathematics Programs. 
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We use our MoSST core dynamics model and geomag- 
netic field at the core-mantle boundary (CMB) continued 
downwarded from surface observations to investigate pos- 
sibilities of geomagnetic data assimilation, so that model 
results and current geomagnetic observations can be used 
to predict geomagnetic secular variation in future. As the 
first attempt, we apply data insertion technique to exam- 
ine evolution of the model solution that is modified by 
geomagnetic input. Our study demonstrate that, with a 
single data insertion, large-scale poloidal magnetic field 
obtained from subsequent numerical simulation evolves 
similarly to the observed geomagnetic variation, regard- 
less the initial choice of the model solution (so long it is a 
well developed numerical solution). The model solution 
diverges on the time scales on the order of 60 years, sim- 
ilar to the time scales of the torsional oscillations in the 
Earth’s core. Our numerical test shows that geomagnetic 
data assimilation is promising with our MoSST model. 


1. Introduction 

Through much of its history, the Earth has possessed 
an internal magnetic field (geomagnetic field) that is be- 
lieved generated and maintained by convective flow in the 
fluid outer core (geodynamo) [Larmor 1919, @]. How- 
ever, it is only for less than a decade, that numerical 
models have been successfully developed to simulate self- 
consistent, fully nonlinear geodynamo processes [Glatz- 
mater and Roberts 1995, Kageyama and Sato 1997, 
Kuang and Bloxham 1997, 6]. These models, though dif- 
ferent in many aspects (e.g. algorithms and physical ap- 
proximations), are able to generate Earth-like magnetic 
field at the CMB. For example, numerical solutions show 
a dominant dipolar field at the CMB, large-scale west- 
ward drift, and occasional field polarity reversals [Konon 
and Roberts 2002, ©]. 

However, geomagnetic and paieomagnetic observe 
tions are not directly utilized in numerical geodynamo 
modeling, mainly because the numerical parameter do- 
mains are far from that for the Earth’s core. This handi- 
caps our understandings on the geodynamo mechanisms, 
thus limiting model improvements and geophysical ap- 
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plications. For example, numerical models can produce 
different solutions that are similar at the CMB, but very 
different deep inside the outer core [Kuang and Bloxham 
1997, Kuang 1999, @]. This is partly caused by dif- 
ferent approximations in the models on torque balances 
on the co-axial (with the Earth’s rotation axis) cylindri- 
cal surfaces across the outer core (the Thylor cylinders). 
Thus surface observations could help identifying appro- 
priate approximations for geodynamo modeling. 

Incorporating observations to numerical modeling can 
facilitate an important application: predicting geomag- 
netic secular variation via data assimilation. This is 
not new, as similar developments occurred in meteo- 
rology and oceanography, where large-scale circulation 
models are used together with past and current observa- 
tions to predict changes in the future. More recently in 
solid Earth research, numerical mantle convection models 
are used together with current observations to hindcast 
historical mantle flow [Bunge et al 2003, Q]. In these 
approaches, appropriate assimilation techniques are se- 
lected to enable us applying the known physics (the mod- 
els) to understand observations (data), and using dis- 
crepancies among model outputs and observational data 
to improve physics knowledges. Similarly in geomag- 
netic data assimilation, observations could be used as 
“time stamps” to modify / constrain the numerical solu- 
tions, such that the modified solutions shall evolve closely 
following the “true” observational trend. 

However, several I im it stl o n - in geomagnetic observa- 
tions could pose serious obstacles to the assimilation. 
First, only the poloidal part Bp of the core field B can be 
observed above the Earth’s surface. The toroidal compo- 
nent Bt is filtered fay a thick, poor electrically conduct- 
ing mantle. Next, Bp is significantly attenuated by the 
crustal magnetic field, leaving only the large scale (for 
degree L < 13 in spherical harmonic expansion) si gnals 
observable above the surface. In addition, data record 
is short and quality decreases back in time: there are 
about 40 years of the highest quality data from global 
satellite measurements [Sa&aJU et al 2002, ©]. Ground 
station and navigation observations provide less accurate 
records over the past centuries [Bloxham and Jackson 
1991, @]. Poorer paleo/archeomagnetic records could ex- 
tend the observed surface field distribution back to more 
than 3000 years [Constable et al 2000, @], Combined the 
observation record is a fraction of the free- decay time 
scale Td (~ 20000 years) in the Earth’s core. Therefore 
the immediate question is whether a modified solution 
with one “time stamp” of geomagnetic observations as- 
similated to a numerical dynamo solution could evolve 
sufficiently close to observations within a reasonable time 
interval. 

Kuang [2000] reported some initial tests on one “time 
stamp” assimilation. His solution suggested that the 
modified solution evolves following much closer to obser- 
vations than purely dynamo simulation. However, no at- 
tempt was made to quantify the tests, in particular error 
development, which is very important for understanding 
the impact of assimilation processes to the dynamics in 
the core, and thus to the evolution of the geomagnetic 
field. 

In this paper we repeat the tests. In particular we 
shall focus on the error development, and the differences 
between the assimilated solutions and the observations. 
The latter shall be used to measure the “improvement” of 
the numerical solutions compared with the observations. 
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2. Assimilation Algorithm 

During the past thirty years, assimilation methods 
have evolved from simple insertion methods, in which 
observation values replace model outputs whenever avail- 
able, to more sophisticated techniques that require de- 
tailed knowledge of error statistics. Since error statistics 
for geodynamo models are unknown, we consider here 
a simple assimilation technique similar to data insertion 
[Berry and Marshall 1989, ©]. Our main purpose is to 
examine the sensitivity of the numerical model to assimi- 
lation process and begin to obtain error information that 
can be used to improve the assimilation scheme. 

This approach can be briefly described as follows. The 
magnetic field B in our model is: 

B = Bt + Bp = Vx(Tlr) + VxVx(Pl,)t 1) 

where 1*. is the radial unit vector, T and P are the 
toroidal and poloidal scalars, respectively. The two 
scalars are expanded in spherical harmonics, 

1<L 

P = Y, «“(»■, t)YT(M) + C.C., (2) 

0 <m<l 

where {Yj m } are fully normalized spherical harmonic 
functions, C.C. denotes the complex conjugate part, and 
(r, 5, <p) define the spherical coordinate. The poloidal 
scalar P is further divided into two parts 

P = Pi + i>2, (3) 

such that 

1<L j 

Pi = ^2 *r y i m + c c. , (4) 

0<m<l 

where Li < L denotes the truncation order in observa- 
tions. 

In the insertion, Pi is first modified at the CMB r = 
fcinb via 

(fcr/tfWim = (&r/&?)ob. = anoTi < u. (5) 

It is then assumed broadcasted “instantly” to the entire 
core via 


(&F7&?)(r< rcn*) = ar- 


Hms, 


1<L\ 

= &?(r) ^2 Q ” y,m + CC - 

0 <m<i 


Next, P 2 is modified as 

PT im = 0P2, 


( 6 ) 


(7) 


( 8 ) 


where the multiplier is determined to conserve the 
poloidal field energy in the outer core 


f [sr * 1 ] 2 dV = f 

J core J coi 


BldV. 


0) 


Other quantities, e.g. the toroidal scalar T and the ve- 
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locity field u, are not changed. 

The modification (7) implies that the insertion is de- 
fined relative to the dipole component &J. The approx- 
imation (6) intends to reduce the transient time far the 
insertion. The conservation of the energy ensures that 
the insertion is neither an energy source nor a sink. 


3. Numerical Results 

In our test here, an arbitrary dynamo solution is se- 
lected to assimilate with the observational data for the 
year 1940. The time evolution of the modified solution 
(the model output) is then used to compare with those 
post 1940 observations derived from the comprehensive 
model [Sabaka et al 2002, Q]. 

In our tests, the parameters are chosen as 

Ro = E = 1.25x10"*, p« = 1, (10) 

we refer the reader to Kuang and Bloxham [1999] for the 
definitions. We also choose L\ = 8 in the insertion, which 
is slightly lower than the order L\ = 13 of the core field 
observed at the Earth’s surface. 

In Figure 1, we plotted several spectral coefficients 
b™/bi of the data insertion test (dashed lines) and of 
the original dynamo simulation (dotted lines) against the 
observations (the solid lines). From the figure we can 
observe that the time evolution of the field coefficients 
have been greatly affected by the insertion: some evolve 
closely to the observations, but some diverge rapidly from 
the observations. However, more coefficients follow the 
observations. In contrast, there is no positive correlation 
between purely dynamo simulation outputs with the ob- 
servations. In fact, & ample error analysis demonstrated 
that the difference between the data insertion outputs 
and the observations are orders of magnitude smaller 
than that between the purely dynamo simulation outputs 
and the observations, as shown in Figure 2. 

Since more coefficients evolve similarly with the obser- 
vations, we expect that the radial component B r of the 
magnetic field at the CMB from the model outputs should 
be similar to those inverted from surface geomagnetic ob- 
servations, as shown in Figure 3. In particular, we can 
observe from the figure that the large scale features, such 
as the south- Atlantic anomaly, are well captured in the 
model outputs, consistent with the coefficients variations. 
However, we should notice that the small scale features, 
such as the magnetic spots in mid Atlantic region (e.g. 
Bermuda area) and in Pacific region are different. While 
some of these features are not consistently reproduced 
within geomagnetic field models, the others indicate that 
continuous assimilation is necessary for better field pre- 
diction. 

The simulation with the one-step data insertion lasts 
approximately 0.1 r* in modeling time (or equivalently 
60 years). After that numerical solutions diverge quickly 
beyond the model resolution limit. Part of the expla- 
nation is that the delicate torque balance on the Taylor 
cylinders is significantly offset by the insertion (6) (see 
Figure 4), resulting in instabilities on the torsional os- 
cillations (with the frequencies of order 0.1 % d) that are 
controlled by the torque balance [Braginsky 1976, 6]. 

4. Discussion 

In this article we described our first attempt on geo- 
magnetic data assimilation. Using a one-step data inser- 
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tion approach (5)-(9), we are able to modify numerical 
outputs such that at the CMB, the time variation of the 
scaled coefficients b^/bi is similar to that of the coef- 
ficients inverted from surface geomagnetic observations. 
Our results suggest that it is possible to assimilate geo- 
magnetic observations into numerical dynamo outputs to 
predict geomagnetic secular variation. 

Our results also reveal several problems that we shall 
work on in the near future for development of a working 
geomagnetic assimilation system. 

First, we shall test the geomagnetic data assimilation 
with a sequence of insertions. This can force the model 
outputs closer to the “true states” with the observed 
geomagnetic secular variation, and provide much better 
knowledge on error development and model sensitivity to 
perturbations on, e.g. the torque balances on the Taylor 
cylinders in the core. Our results show that a strong per- 
turbation on the torque balance could quickly bring the 
model outputs close to the observations, but destabilize 
the simulation in a short period. A weaker perturbation 
may reduce instability problem, however slow down the 
assimilation. Understanding the error development and 
model sensitivity will require farther studies using ob- 
servation similation experiments (OSE’s) and Newtonian 
nudging [Davies and Turner 1977, ©] that could eliminate 
some problems, e.g. non-physical oscillations, inherited 
from the direct insertion technique. OSE’s involve the 
use of two models (or one model with different parame- 
ter values), one of which is treated as the "truth” and the 
other is the model. Observations taken from the former 
are then assimilated into the latter. 

Another problem to be addressed soon is the appropri- 
ate scale between the modeling time and the observation 
time. The best scaling rules depend on the torsional os- 
cillations which have the frequencies ~ i^ 1/2 and decide 
short-time secular variations [Braginsky 1976, 6]. How- 
ever, with the parameters (10), the torsional oscillations 
are not separable from other dynamo waves controlled 
by the leading order magnetostrophic balance in the core 
[Kuang and Bloxham 1999, Q]. Therefore, we need to 
identify a proper parameter domain in which the tor- 
sional oscillations are well separated from other dynamo 
waves, and computation demand is still manageable. 
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Figure 1. Time variation of the spectral coefficients 
by 1 /bi at the CMB from observations (solid lines), data 
assimilation (dashed lines) and numerical dynamo mod- 
eling (dotted lines). The blue and red lines are the real 
and imaginary parts of 6J”, respectively. The green lines 
are for 63. 
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Figure 2. Snapshots of the error analysis of the spectral 
coefficients 6™ from observations and the model outputs. 
The top panel are the errors in the magnitude, and the 
bottom panel are the erros in the phase. The left panel 
are the results without insertion, and the right are the 
results with the insertion (assimilation). 
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Figure 3. Snapshots of the radial component B r at 
the CMB from surface observations (left panel) and from 
assimilation (right panel). 
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Figure 4. The distributions of the axial torque T* = 
/ E ( J X B)+ dS (J is the current density and the subscript 
<f> denotes the zonal component) on the Thylor cylinders 
£ across the the outer core. The dotted line is 50r, for 
the dynamo solution, and the solid line is T z modified by 
the data insertion. 




